Data

This study uses landslide, fire, and precipitation data to study the antecedent precipitation of post-wildfire landslides. The impact of wildfire on rainfall-triggered landslide susceptibility in a variety of global regions is evaluated through a proxy of normalized precipitation.

Data Filepaths

Use this section to modify file paths if running the notebook on a different computer.

# Landslide data from the NASA Global Landslide Catalog
slide.path <- '../regions_data/glc/GLC20180821.csv'

# Fire data from the MODIS Burned Area dataset
modis.paths <- list.files(path='../regions_data/processed', 
                          pattern='slide_modisburndate_[ns][ew].csv',
                          full.names = T)

# Precipitation data from CHIRPS
precip.chirps.ww.path <-  file.path("../regions_data/processed", 
                                    "glc_precip_global.csv")
precip.chirps.ww.monthly.path <- file.path("../regions_data/", "processed",
                                   "glc_precip_global_monthly.csv")

# Additional preciptiation data from Daymet for a limited region
precip.daymet.ca.path <- file.path("../regions_data/processed", 
                                   "glc_precip_daymet.csv")

Data inclusion

Include landslides with the following parameters: * Rainfall triggered

rain.triggers <- c('rain', 'downpour', 'flooding', 'continuous_rain')
  • Occurred within the time frame of the fire and precipitation datasets
start.date <- lubridate::ymd('2000-01-01')
end.date <- lubridate::ymd('2019-01-01')
  • Occurred within the spatial extend of the fire and precipitation datasets
min.latitude <- -50
max.latitude <- 50
  • Location accuracy on par with precipitation resolution or more accurate
min.location.accuracy <- '10km'
  • MODIS Burned Area data should include an additional 3 years
modis.start.date <- lubridate::ymd('20040101')
modis.end.date <- lubridate::ymd('20190101')
  • Wilcox tests between burned and unburned locations have the following parameters:
wilcox.alternative <- 'less' # n - burned, m - unburned
wilcox.alpha <- 0.05

Study definitions

  • Study location parameters are applied using a set of functions to identify if data are included in the study.
# Identify if a location is in the global limits
is.in.global.study <- function (x, y) { between(y, min.latitude, max.latitude) }

# Identify if a location is in the CA/NV pilot study
get.data('state.boundaries', tigris::states)
boundaries <- subset(state.boundaries, STUSPS %in% c('CA', 'NV'))
boundaries.df <- fortify(boundaries)
is.in.ca.study <- function (x, y) {
  bbox <- sf::st_bbox(subset(state.boundaries, STUSPS == 'CA'))
  between(x, bbox$xmin, bbox$xmax) &
    between(y, bbox$ymin, bbox$ymax)
}
  • Seasons are split in the Northern Hemisphere by:

    Winter
    December-January-February
    Spring
    March-April-May
    Summer
    June-July-August
    Fall
    September-October-November
  • In the southern hemisphere, the split is shifted instead:

    Winter
    June-July-August
    Spring
    September-October-November
    Summer
    December-January-February
    Fall
    March-April-May
  • A set of functions converts months and days-of-the-year to seasons.

season.breaks <- c('spring'=59, 'summer'=151, 'fall'=243, 'winter'=334)
yday.to.season <- function(yday, latitude) {
    if (latitude > 0) {
      case_when(
        yday > season.breaks['winter'] ~ 'Winter',
        yday > season.breaks['fall'] ~ 'Fall',
        yday > season.breaks['summer'] ~ 'Summer',
        yday > season.breaks['spring'] ~ 'Spring',
        TRUE ~ 'Winter')
    } else {
      case_when(
        yday > season.breaks['winter'] ~ 'Summer',
        yday > season.breaks['fall'] ~ 'Spring',
        yday > season.breaks['summer'] ~ 'Winter',
        yday > season.breaks['spring'] ~ 'Fall',
        TRUE ~ 'Summer')
    }
}

month.to.season <- function (month, north) {
  if (north) {
    case_when(
      month >= 12 ~ 'Winter',
      month >= 9 ~ 'Fall',
      month >= 6 ~ 'Summer',
      month >= 3 ~ 'Spring',
      TRUE ~ 'Winter')
  } else {
    case_when(
      month >= 12 ~ 'Summer',
      month >= 9 ~ 'Spring',
      month >= 6 ~ 'Winter',
      month >= 3 ~ 'Fall',
      TRUE ~ 'Summer')
  }
}
  • Consistent plot colors and themes
# Region colors
region.plot.order <- c('All landslides',
                       'California',
                       'Intermountain West',
                       'Pacific Northwest',
                       'Central America',
                       'Himalayas',
                       'Southeast Asia')
scale_region_values <- c('All landslides' = 'grey50',
                         setNames(RColorBrewer::brewer.pal(6, 'Dark2'), 
                                  region.plot.order[2:7]))
scale_region_labels <- c('All landslides' = 'None',
                         setNames(region.plot.order[2:7], 
                                  region.plot.order[2:7]))
scale_color_region <- scale_color_manual('Region', 
                                         values = scale_region_values,
                                         labels = scale_region_labels)
scale_fill_region <- scale_fill_manual('Region', 
                                       values = scale_region_values,
                                       labels = scale_region_labels)

scale_color_precip <- scale_color_distiller('Precipitation percentile', 
                                            palette='BrBG', direction = 1,
                                            guide = guide_colorbar(barwidth = unit(0.5, 'npc')))
scale_fill_precip <- scale_fill_distiller('Precipitation percentile', 
                                          palette='BrBG', direction = 1,
                                          guide = guide_colorbar(barwidth = unit(0.5, 'npc')))

# Fire colors
burned.labels <- c('burned'='Burned', 'unburned'='Unburned')
burned.colors <- c('burned' = '#d95f02', unburned = '#7570b3')
scale_color_burned <- scale_color_manual('', 
                                         values = burned.colors,
                                         labels = burned.labels)
scale_fill_burned <- scale_fill_manual('', 
                                       values = burned.colors,
                                       labels = burned.labels) 
# Season
scale_x_season <- scale_x_continuous(
    breaks = season.breaks,
    labels = c('Spring', 'Summer', 'Fall', 'Winter'))

col.ocean <- 'lightskyblue1'
theme_map <- list(
  theme_bw(),
  coord_quickmap(),
  labs(x='', y='')
)

theme_global_map <- list(theme_map,
                         borders(fill = 'grey90', size = 0.1),
                         theme(panel.background = element_rect(fill=col.ocean), 
                               panel.grid = element_blank(),
                               axis.text = element_blank(),
                               axis.ticks = element_blank(),
                               plot.margin = unit(c(0,0,0,0), "null")),
                         coord_quickmap(ylim = c(-60, 75), xlim = c(-180, 180))
)
theme_panel <- theme(panel.background = element_rect(fill=col.ocean), 
                     panel.grid = element_blank(),
                     axis.text = element_blank(),
                     axis.ticks = element_blank(),
                     plot.margin = unit(c(1,0,-1,0), "lines"),
                     legend.position = 'null')
theme_panel_legend <- theme(
  legend.box.margin = margin(-1, 0, 1, 0),
  legend.box.just = c(0, 0),
  #legend.background = element_rect(color='black'), #debug
  legend.position = 'bottom'
)
capitalize <- function(string) {
  substr(string, 1, 1) <- toupper(substr(string, 1, 1))
  string
}

Landslide data from the NASA Global Landslide Catalog

load.slides <- function (slide.path, 
                         start.date, end.date, 
                         is.in.global.study, 
                         rain.triggers, 
                         min.location.accuracy) {
  # Define factor levels
  ls.sizes <- c("small", "medium", "large", "very_large", "catastrophic")
  ls.accuracy <- c("exact", "1km", "5km", "10km", 
                   "25km", "50km", "100km", "250km")
  
  # Read data file and remove duplicates
  slide.raw <- readr::read_csv(slide.path,
    guess_max = 100000,
    na = c("", "unknown", "--"),
    col_types = cols(
      landslide_size = readr::col_factor(ls.sizes, ordered=TRUE),
      location_accuracy = readr::col_factor(ls.accuracy, ordered=TRUE),
      landslide_trigger = readr::col_factor(NULL),
      landslide_category = readr::col_factor(NULL),
      landslide_setting = readr::col_factor(NULL),
      country_code = readr::col_factor(NULL)
    )
  ) %>% 
    dplyr::mutate(event_date = as_date(ymd_hm(event_date)),
                  submitted_date = as_date(ymd_hm(submitted_date)),
                  last_edited_date = as_date(ymd_hm(last_edited_date)))
  slide.unique <- slide.raw  %>%
    dplyr::group_by(event_date, latitude, longitude) %>%
    dplyr::mutate(i = row_number()) %>%
    dplyr::filter(i == 1) %>%
    dplyr::select(-i) %>%
    dplyr::ungroup()
  
  # Filter landslides to meet study parameters
  slide.unique %>%
    dplyr::filter(event_date >= start.date,
                  event_date <= end.date,
                  is.in.global.study(longitude, latitude),
                  landslide_trigger %in% rain.triggers,
                  location_accuracy <= min.location.accuracy) %>%
    droplevels() %>%
    dplyr::mutate(event_yday = yday(event_date),
                  debris_flow = landslide_category %in% c('mudslide', 'debris flow'),
                  season = factor(map2_chr(event_yday, latitude, 
                                           ~yday.to.season(.x, .y)),
                                  levels = c('Spring', 'Summer', 'Fall', 'Winter'),
                                  ordered = T))
}

get.data('slide.all.precip', overwrite = F,
         load.slides,
         slide.path, start.date, end.date, 
         is.in.global.study, rain.triggers,
         min.location.accuracy)

slide.all.precip %>%
  dplyr::select(landslide_trigger, location_accuracy, latitude, event_date) %>%
  summary()
##        landslide_trigger location_accuracy    latitude     
##  rain           :1823    exact: 466        Min.   :-46.77  
##  continuous_rain: 588    1km  :1707        1st Qu.: 10.49  
##  downpour       :3254    5km  :2437        Median : 29.91  
##  flooding       :  48    10km :1103        Mean   : 23.52  
##                                            3rd Qu.: 39.11  
##                                            Max.   : 49.98  
##    event_date        
##  Min.   :2000-10-13  
##  1st Qu.:2010-11-30  
##  Median :2013-06-24  
##  Mean   :2013-05-17  
##  3rd Qu.:2016-03-11  
##  Max.   :2018-08-07

Precipitation data from CHIRPS and Daymet

Global CHIRPS data

load.precip.all <- function (precip.chirps.ww.path) {
    precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
                                     col_types = cols(X1 = col_datetime(format = "%Y-%m-%d")))
    precip.chirps.ww.raw %>%
      dplyr::select(date=X1, OBJECTID, pctl=precip_pctl_7) %>%
      dplyr::filter(OBJECTID %in% slide$OBJECTID) %>%
      dplyr::mutate(window=7) %>%
      dplyr::filter(!is.na(pctl)) # Remove NA values from beginning of data range
}

get.data('precip.chirps.ww.all', precip.chirps.ww.path)
precip.chirps.ww.all

Detailed global CHIRPS data close to the landslide

date.filter.precip <- function (raw, slide=slide.all.precip, pre=90, post=30) {
  raw %>%
    dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date), 
                      by='OBJECTID') %>% 
    dplyr::select(OBJECTID, dplyr::starts_with('precip_'), event_date, date) %>%
    dplyr::mutate(date = lubridate::as_date(date),
           days_before = as.integer(event_date - date)) %>%
    dplyr::filter(days_before < pre, days_before > -post)
}

tidy.precip <- function (precip) {
  precip %>% 
    dplyr::select(OBJECTID, starts_with('precip_'), event_date, date) %>%
    gather(key='variable', value='amount', starts_with('precip_')) %>%
    separate(variable, into=c(NA, 'variable', 'window')) %>% 
    dplyr::filter(variable=='pctl') %>%
    dplyr::select(-variable, pctl=amount) %>%
    dplyr::mutate(days_before = as.integer(event_date - date),
           yday = yday(event_date - days(days_before)),
           window = as.numeric(window)) %>%
    dplyr::select(-event_date)
}

load.event.precip <- function (precip.chirps.ww.path) {
    precip.chirps.ww.raw <- read_csv(precip.chirps.ww.path,
                                     col_types  = readr::cols(X1 = col_datetime(format = "%Y-%m-%d")))
    precip.chirps.ww.raw.filt <- precip.chirps.ww.raw %>%
      dplyr::rename(date = X1) %>%
      dplyr::filter(year(date) > 2003)
    precip.chirps.ww.raw.filt2 <- precip.chirps.ww.raw.filt %>%
      date.filter.precip 
    precip.chirps.ww.raw.filt2 %>% 
      tidy.precip
}

get.data('precip.chirps.ww.event', precip.chirps.ww.path)
precip.chirps.ww.event

Import Monthly CHIRPS data

load.precip.monthly <- function (monthly.path, slide) {
  read_csv(precip.chirps.ww.monthly.path,
           col_types = cols(X1 = col_date(format = '%Y-%m-%d'))) %>%
    dplyr::transmute(OBJECTID, 
                     year = year(X1), month = month(X1), 
                     precip.mm=precip) %>%
    dplyr::inner_join(slide %>% select(OBJECTID, latitude), by='OBJECTID') %>%
    dplyr::mutate(north = latitude > 0) %>%
    dplyr::group_by(north) %>%
    dplyr::mutate(season = factor(month.to.season(month, north),
                                  levels = c('Spring', 'Summer', 'Fall', 'Winter'),
                                  ordered = T)) %>%
    ungroup()
}
 
get.data('precip.chirps.ww.monthly', overwrite = F,
         load.precip.monthly,
         precip.chirps.ww.monthly.path,
         slide.all.precip)

precip.chirps.ww.seasonal.median <- precip.chirps.ww.monthly %>%
  group_by(OBJECTID, year, season) %>%
  summarize(precip.mm = sum(precip.mm)) %>%
  ungroup() %>%
  summarize(median = median(precip.mm)) %>%
  pull(median)

precip.chirps.ww.monthly

Filter out landslides with zero precipitation

precip.chirps.ww.event %>%
  # Filter out sites with no precipitation data at all
  filter(days_before==-1, window==7, !is.na(pctl)) %>%
  # Filter out sites already excluded by landslide criteria
  inner_join(slide.all.precip, by = 'OBJECTID') %>%
  # Count numbers of zero and non-zero precipitation
  group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
  tally() %>%
  spread(key = precip, value = n) %>%
  mutate(total = non.zero.precip + zero.precip,
         percent.zero.precip = zero.precip / total)
filter.slide.precip <- function (slide.all.precip, precip.event) {
  slide.all.precip %>%
    dplyr::inner_join(precip.event %>% 
                        filter(days_before==-1, window==7, !is.na(pctl)) %>%
                        dplyr::filter(days_before==-1, window==7, pctl > 0)) %>%
    dplyr::select(OBJECTID, latitude, longitude, event_date, event_yday, season,
                  location_accuracy, 
                  landslide_category, landslide_trigger, landslide_size,
                  landslide_setting, debris_flow)
}
get.data('slide',  overwrite = F,
         filter.slide.precip,
         slide.all.precip,
         precip.chirps.ww.event)
slide

Filter California Landslides

filter.slide.ca <- function (slide.all.precip) {
  slide.all.precip %>%
    filter(is.in.ca.study(longitude, latitude)) %>%
    mutate(region = as.factor(case_when(
             latitude < 35.0 ~ 'South',
             longitude < -121.25 ~ 'North',
             TRUE ~ 'East')))
}

get.data('slide.ca',  overwrite = F,
         filter.slide.ca, 
         slide.all.precip)

slide %>% summary
##     OBJECTID         latitude        longitude          event_date        
##  Min.   :388907   Min.   :-46.77   Min.   :-179.865   Min.   :2007-01-05  
##  1st Qu.:391763   1st Qu.: 10.33   1st Qu.: -96.159   1st Qu.:2010-11-30  
##  Median :394578   Median : 28.55   Median :  -8.624   Median :2013-06-16  
##  Mean   :394614   Mean   : 23.09   Mean   :  -2.247   Mean   :2013-05-15  
##  3rd Qu.:397488   3rd Qu.: 38.95   3rd Qu.:  92.212   3rd Qu.:2016-03-10  
##  Max.   :400281   Max.   : 49.98   Max.   : 179.402   Max.   :2018-08-07  
##                                                                           
##    event_yday       season     location_accuracy   landslide_category
##  Min.   :  1.0   Spring:1310   exact: 423        landslide  :3492    
##  1st Qu.: 91.0   Summer:1798   1km  :1573        mudslide   :1326    
##  Median :172.0   Fall  :1104   5km  :2259        rock_fall  : 226    
##  Mean   :171.6   Winter:1101   10km :1058        complex    : 103    
##  3rd Qu.:244.0                                   debris_flow:  98    
##  Max.   :366.0                                   (Other)    :  64    
##                                                  NA's       :   4    
##        landslide_trigger      landslide_size     landslide_setting
##  rain           :1656    small       :1735   above_road   :1257   
##  continuous_rain: 561    medium      :3160   natural_slope: 265   
##  downpour       :3054    large       : 326   urban        : 178   
##  flooding       :  42    very_large  :  36   below_road   : 121   
##                          catastrophic:   3   above_river  :  65   
##                          NA's        :  53   (Other)      : 187   
##                                              NA's         :3240   
##  debris_flow    
##  Mode :logical  
##  FALSE:3987     
##  TRUE :1326     
##                 
##                 
##                 
## 

Filter California CHIRPS data

filter.precip.ca <- function (precip.chirps.ww.event, slide.ca) {
  precip.chirps.ca <- precip.chirps.ww.event%>%
    dplyr::filter(OBJECTID %in% slide.ca$OBJECTID)
}

get.data('precip.chirps.ca.event', filter.precip.ca,  overwrite = F,
         precip.chirps.ww.event, slide.ca)
get.data('precip.chirps.ca.all', filter.precip.ca, overwrite = F,
         precip.chirps.ww.all, slide.ca)

precip.chirps.ca.event %>%
  # Filter out sites with no precipitation data at all
  filter(days_before==-1, window==7, !is.na(pctl)) %>%
  # Filter out sites already excluded by landslide criteria
  inner_join(slide.all.precip, by = 'OBJECTID') %>%
  # Count numbers of zero and non-zero precipitation
  group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
  tally() %>%
  spread(key = precip, value = n) %>%
  mutate(total = non.zero.precip + zero.precip,
         percent.zero.precip = zero.precip / total * 100)

Import Daymet data for California

load.precip.daymet.ca <- function (precip.daymet.ca.path) {
  read_csv(precip.daymet.ca.path, 
           col_types = cols(X1 = col_datetime(format = "%Y-%m-%d %H:%M:%S"))) %>% 
  rename(date=X1) %>%
  date.filter.precip %>%
  tidy.precip
}

get.data('precip.daymet.ca.all', overwrite = F,
         load.precip.daymet.ca, 
         precip.daymet.ca.path)

precip.daymet.ca.all

Fire data computed from MODIS Burned Area

load.modis <- function (modis.paths, burn.date.min, burn.date.max) {
  modis.date.raw <- dplyr::bind_rows(map(modis.paths, read_csv))
  burn.record.length <- time_length(burn.date.max - burn.date.min, 'years')
  modis.date.raw %>% 
    dplyr::right_join(slide %>% dplyr::select(OBJECTID, event_date, latitude), 
                      by='OBJECTID') %>%
    # Convert MODIS Burned Area value to date
    dplyr::mutate(burn_date = make_date(year=floor(month / 1000)) + burn_date) %>%
    dplyr::select(-month) %>%
    # Separate out indivicual fires
    dplyr::group_by(OBJECTID) %>%
    dplyr::mutate(gap = as.integer(burn_date - lag(burn_date, default=min(burn_date))),
                  new.fire = gap > 30,
                  fire.id = cumsum(new.fire)) %>%
    drop_na(fire.id) %>%
    # Calculate burn fraction, fire ending date, and season for each fire
    dplyr::group_by(OBJECTID, event_date, fire.id) %>%
    summarize(fraction = sum(fraction), 
              burn_date = max(burn_date),
              burn_season = yday.to.season(yday(burn_date), latitude)) %>%
    # Calculate time relative to landslide
    dplyr::mutate(burn_days_before = event_date - burn_date,
                  burn_years_before = time_length(burn_days_before, 'years')) %>%
    # Compute fire information
    dplyr::group_by(OBJECTID) %>%
    dplyr::mutate(pre.slide = burn_date==max(burn_date[burn_days_before>=0]),
                  n = n(),
                  period = burn.record.length / n,
                  start_days_before = min(event_date) - burn.date.min,
                  start_years_before = time_length(start_days_before, 'years'),
                  end_days_before = max(event_date) - burn.date.max,
                  end_years_before = time_length(end_days_before, 'years')) %>%
    dplyr::select(-event_date) %>%
    dplyr::ungroup()
}

get.data('modis.date', overwrite = F,
         load.modis,
         modis.paths, modis.start.date, modis.end.date)
modis.date
classify.burned <- function (modis.date, slide) {
  modis.date  %>%
    dplyr::filter(burn_years_before < 3, pre.slide) %>%
    dplyr::right_join(slide, by='OBJECTID') %>%
    dplyr::group_by(OBJECTID) %>%
    summarize(burn.cumulative = sum(fraction)) %>%
    dplyr::mutate(burn.cumulative = ifelse(is.na(burn.cumulative), 
                                           0, 
                                           burn.cumulative),
                  burned = ifelse(burn.cumulative > 0, 'burned', 'unburned')) %>%
    dplyr::select(OBJECTID, burned, burn.cumulative) %>% 
    dplyr::ungroup()
}

get.data('modis', overwrite = F,
         classify.burned, 
         modis.date,
         slide)

modis %>%
  group_by(burned) %>%
  tally() %>%
  mutate(percent = n / nrow(slide) * 100)

Define regions using AGNES clustering

cluster.global <- function (slide, agclusts, k=30) {
  # Select clusters with more than 100 landslides, and label
  slide %>%
    dplyr::select(OBJECTID, latitude, longitude) %>%
    dplyr::mutate(region = agclusts %>% cutree(k)) %>%
    dplyr::group_by(region) %>%
    mutate(region = as_factor(region)) %>%
    dplyr::filter(n() > 100)  %>%
    dplyr::ungroup() %>%
    dplyr::mutate(region = fct_recode(as_factor(region),
                              'Central Europe' = '1', 
                              'Malaysia' = '2', 
                              'Western US/Canada' = '3',
                              'East Brazil' = '4',
                              'Southeast Asia' = '5',
                              'South India' = '6',
                              'Eastern US' = '8',
                              'Himalayas' = '9',
                              'Southeast Asia' = '10',
                              'Caribbean/Venezuela' = '11', 
                              'Central America' = '14'))
}

# Cluster landslides based on location using AGNES algorithm
get.data('agclusts', 
         cluster::agnes, 
         slide %>%
           dplyr::select(latitude, longitude))


get.data('clusters.global',  overwrite = F,
         cluster.global, 
         slide.all.precip, agclusts)

clusters.global %>%
    ggplot(aes(longitude, latitude)) +
    theme_global_map +
    geom_point(aes(color = region), size=1) +
    scale_color_brewer('Region', palette='Paired') +
    guides(color=guide_legend(override.aes = list(size=5)))
split.regions.us <- function (agclusts, slide, k=6) {
  # Crop cluster tree, select clusters with more than 100 landslides, and label
  slide %>%
    dplyr::select(OBJECTID, latitude, longitude) %>%
    dplyr::mutate(region = agclusts %>% cutree(k)) %>%
    mutate(region = as_factor(region)) %>%
    dplyr::mutate(region = fct_recode(as_factor(region),
                                      'Intermountain West' = '3',
                                      'Intermountain West' = '4',
                                      'Intermountain West' = '5',
                                      'Pacific Northwest' = '1',
                                      'California' = '2')) %>%
    dplyr::group_by(region) %>%
    dplyr::filter(n() > 100)  %>%
    dplyr::ungroup() %>%
    droplevels()
}

get.data('agclusts.us',  overwrite = F,
         function (coord) cluster::agnes(coord, stand=T), 
         clusters.global %>%
           dplyr::filter(region=='Western US/Canada') %>%
           dplyr::select(latitude, longitude))

get.data('clusters.us',  overwrite = F,
         split.regions.us,
         agclusts.us, 
         clusters.global %>%
           dplyr::filter(region=='Western US/Canada'))

clusters.us %>%
  ggplot(aes(longitude, latitude)) +
  geom_point(aes(color = region), size=1) +
  scale_color_brewer(palette='Set2') +
  guides(color=guide_legend(override.aes = list(size=5))) +
  theme_map
merge.regions <- function (clusters.global, clusters.us) {
  # Merge US and global regions
  clusters.global %>%
    left_join(clusters.us %>% dplyr::select(OBJECTID, region), 
              by='OBJECTID') %>%
    dplyr::mutate(region = as_factor(ifelse(!is.na(region.y), 
                                            as.character(region.y), 
                                            as.character(region.x)))) %>%
    dplyr::select(-starts_with('region.')) %>%
    filter(region != 'Western US/Canada') %>%
    droplevels()
}

get.data('regions.all', overwrite = F,
         merge.regions, 
         clusters.global, clusters.us)
get.regions.info <- function (regions.all, slide) {
  regions.all %>%
    bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>%
    dplyr::inner_join(modis, by = 'OBJECTID') %>%
    group_by(region) %>%
    dplyr::group_by(region, burned) %>%
    dplyr::summarize(n = n()) %>%
    tidyr::spread(key = burned, value = n) %>%
    dplyr::mutate(n = burned + unburned,
                  percent.burned = burned / n * 100) %>%
    ungroup
}

get.data('regions.info', overwrite = F,
         get.regions.info,
         regions.all, slide)

regions.info
filter.and.merge.regions <- function (regions.all, regions.info, slide) {
  # Combine nearby regions
  regions.all  %>%
    filter(OBJECTID %in% slide$OBJECTID) %>%
    bind_rows(slide %>% dplyr::mutate(region = 'All landslides')) %>% 
    dplyr::left_join(regions.info) %>% 
    mutate(region = fct_recode(region,
                               'Central America' = 'Caribbean/Venezuela')) %>%
    dplyr::filter(region != 'East Brazil', 
                  percent.burned > 4) %>%
    droplevels() %>%
    dplyr::mutate(region = factor(region, region.plot.order, ordered = T)) %>%
    dplyr::select(OBJECTID, region)
}

get.data('slide.regions',  overwrite = F,
         filter.and.merge.regions, 
         regions.all, 
         regions.info, 
         slide)

regions.tally <- slide.regions %>%
  group_by(region) %>%
  tally

regions.tally

Exploration

Map

A panelled maps showing landslide regions (a and b), whether or not each landslide took place on a burned site (c and d), and the precipitation percentile of the landslide-triggering storm (e and f). Inset plots show, for each region, the distribution of burned fractions of the landslide sites (a and b), location accuracy of landslides split by fire history (c and d), and the seasonal total precipitation climatology (e and f). ### Regions map

gg.region <- slide.regions %>%
  left_join(slide, by='OBJECTID') %>%
  mutate(region = fct_relevel(region, 'All landslides', after = Inf)) %>%
  arrange(desc(region)) %>%
  ggplot(aes(longitude, latitude)) +
  theme_global_map +
  geom_point(aes(color = region), size=0.02) + 
  labs(color = '') +
  scale_color_region +
  guides(color = guide_legend(override.aes = list(size = 2))) +
  theme_panel_legend

gg.region

Burned sites map

gg.burned <- slide %>%
  dplyr::left_join(modis) %>%
  dplyr::arrange(desc(burned)) %>%
  ggplot() +
  theme_global_map +
  #geom_point(data = slide.regions %>%
  #             dplyr::filter(region !='All landslides') %>%
  #             dplyr::left_join(slide, by='OBJECTID'),
  #           aes(x = longitude, y = latitude, fill = region),
  #           size = 2, alpha = .01, shape = 21, stroke = NA) +
  geom_point(aes(x=longitude, y=latitude, color=burned),
             size=.02) +
  scale_color_burned +
  #scale_fill_region +
  guides(color = guide_legend(override.aes = list(size = 2),
                              title = NULL),
         fill = guide_none()) +
  theme_panel_legend

gg.burned

Precipitation map

gg.precip <- slide %>%
  left_join(precip.chirps.ww.event %>% 
              filter(days_before==-1, window==7), 
            by='OBJECTID') %>%
  filter(pctl > 0) %>%
  arrange(desc(pctl)) %>%
  ggplot() +
  theme_global_map +
  geom_point(aes(x=longitude, y=latitude, color=pctl),
             size=.02) +
  scale_color_precip +
  theme_panel_legend

gg.precip

Burned fraction

min.burn.cumulative <- modis %>%
  filter(burn.cumulative > 0) %>%
  pull(burn.cumulative) %>%
  min
max.burn.cumulative <- modis %>%
  pull(burn.cumulative) %>%
  max
plot.burned.fraction <- function (reg, slide.regions, modis) {
  slide.regions %>%
    filter(region==reg) %>%
    left_join(modis, by = 'OBJECTID') %>%
    filter(burn.cumulative > 0) %>%
    ggplot() +
    geom_density(aes(x = burn.cumulative), fill = 'black') +
    scale_x_log10() +
    labs(x = 'Burned fraction', y = '') +
    scale_color_region +
    scale_fill_region +
    coord_cartesian(expand=F, xlim = c(min.burn.cumulative, max.burn.cumulative)) +
    theme_bw() +
    theme(legend.position = 'none',
          axis.text.x = element_text(angle = 90, vjust = 0.5, size = 4),
          axis.title = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          plot.background = element_blank(),
          panel.grid = element_blank(),
          plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
    ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}

gg.burned.fraction <- region.plot.order %>%
  set_names %>%
  map(plot.burned.fraction, slide.regions, modis)

cowplot::plot_grid(plotlist = gg.burned.fraction, nrow = 1)

Location accuracy by region/fire history

plot.accuracy <- function (reg, slide.regions, slide, modis) {
  slide.regions %>%
    left_join(slide, by='OBJECTID') %>%
    left_join(modis, by='OBJECTID') %>%
    filter(region==reg) %>%
    mutate(location_accuracy = fct_recode(location_accuracy, '0km' = 'exact')) %>%
    drop_na(location_accuracy) %>%
    ggplot() +
    geom_bar(aes(x = location_accuracy, fill=burned,
                 group = interaction(region, burned),
                 y = ..prop..), 
             position=position_dodge2(preserve = 'single')) +
    scale_fill_burned +
    coord_cartesian(expand=F, ylim = c(0, 0.6)) +
    theme_bw() +
    labs(y = '', x = '', fill = 'Location\naccuracy') +
    theme(legend.position = 'none',
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, 
                                     size = 4),
          axis.title = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          plot.background = element_blank(),
          panel.grid = element_blank(),
          plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
    ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}

gg.accuracy <- region.plot.order %>%
  set_names %>%
  map(plot.accuracy, slide.regions, slide, modis)

cowplot::plot_grid(plotlist = gg.accuracy, nrow = 1)

Climatology by region

plot.climatology <- function (reg, precip.monthly, slide.regions, med) {
  precip.seasonal <- precip.monthly %>%
    group_by(OBJECTID, year, season) %>%
    summarize(precip.mm = sum(precip.mm)) 
  
  ymin <- 1
  ymax <- max(precip.seasonal$precip.mm, na.rm=T)
  
  precip.seasonal %>%
    dplyr::right_join(slide.regions %>% filter(region==reg)) %>%
    drop_na(season) %>%
    ggplot()  +
    geom_violin(aes(y = precip.mm, x = season, group = season),
                fill='black', color=NA, adjust = 2)  +
    geom_hline(aes(yintercept = med), color = 'grey', size = .3)+
    scale_y_log10(breaks = c(1, 10, 100, 1000)) +
    scale_x_discrete(labels = c('Spring' = 'MAM', 'Summer' = 'JJA',
                                'Fall' = 'SON', 'Winter' = 'DJF'),
                     expand = ggplot2::expand_scale(add = c(1, 0))) +
    coord_cartesian(expand=F, ylim = c(ymin, ymax)) +
    theme_bw() +
    theme(legend.position = 'none',
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 4, 
                                     margin=margin(0,0,0,0)),
          axis.text.y = element_text(size = 4, hjust = 1, vjust = 0.5),
          axis.title = element_blank(),
          plot.background = element_blank(),
          panel.grid = element_blank(),
          plot.title = element_text(size = 5, margin=margin(0,0,0,0))) +
    ggtitle(str_c(reg, '\nn=', pull(filter(regions.tally, region==reg))))
}

gg.climatology <- region.plot.order %>%
  purrr::set_names() %>%
  purrr::map(plot.climatology, 
             precip.chirps.ww.monthly, 
             slide.regions, 
             precip.chirps.ww.seasonal.median)
 
cowplot::plot_grid(plotlist = gg.climatology, nrow = 1)

Full map figure

coord_america <- coord_equal(xlim = c(-140, -50), ylim = c(-8, 53), expand = F)
coord_asia <-  coord_equal(xlim = c(50, 140), ylim = c(-12, 49), expand = F) 

# Clip maps
gg.america.region <- gg.region + theme_panel + coord_america
gg.asia.region <- gg.region + theme_panel + coord_asia
gg.america.burned <- gg.burned + theme_panel + coord_america
gg.asia.burned <- gg.burned + theme_panel + coord_asia
gg.america.precip <- gg.precip + theme_panel + coord_america
gg.asia.precip <- gg.precip + theme_panel + coord_asia

# Pull legends
legend.regions <- cowplot::get_legend(gg.region)
legend.burned <- cowplot::get_legend(gg.burned)
legend.precip <- cowplot::get_legend(gg.precip)

# Region midpoints for arrows
region.centroids <- slide.regions %>%
  left_join(slide) %>%
  group_by(region) %>%
  summarize(y = mean(latitude), x = mean(longitude))
# Add insets
inset.w <- .2
inset.h <- .41
inset.arrow <- grid::arrow(length = unit(0.03, "npc"), type = 'open')
# Regions
gg.america.region.inset <- cowplot::ggdraw(gg.america.region) +
  # Arrows
  cowplot::draw_grob(grid::grid.lines(x = c(.1, .25), y = c(.8, .8), 
                                      arrow=inset.arrow)) +
  cowplot::draw_grob(grid::grid.lines(x = c(.1, .275), y = c(.4, .65), 
                                      arrow=inset.arrow)) +
  cowplot::draw_grob(grid::grid.lines(x = c(.4, .4), y = c(.375, .575), 
                                      arrow=inset.arrow)) +
  cowplot::draw_grob(grid::grid.lines(x = c(.825, .675), y = c(.675, .3), 
                                      arrow=inset.arrow)) +
  # Inset plots
  cowplot::draw_plot(gg.accuracy$`Pacific Northwest`, 
                     x = 0.01, y = 1, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.accuracy$California, 
                     x = 0.01, y = .6, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.accuracy$`Intermountain West`, 
                     x = .275, y = .4, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.accuracy$`Central America`,  
                     x = .775, y = .8, 
                     width = inset.w, height = inset.h, vjust = 1)

gg.asia.region.inset <- cowplot::ggdraw(gg.asia.region) +
  #Arrows
  cowplot::draw_grob(grid::grid.lines(x = c(.3, .4), y = c(.3, .6), 
                                      arrow=inset.arrow)) +
  cowplot::draw_grob(grid::grid.lines(x = c(.825, .8), y = c(.55, .4), 
                                      arrow=inset.arrow)) +
  # Inset plots
  cowplot::draw_plot(gg.accuracy$Himalayas,
                     x = .15, y = 0.4, 
                     width = inset.w, height= inset.h, vjust = 1) +
  cowplot::draw_plot(gg.accuracy$`Southeast Asia`,
                     x = .775, y = .935, 
                     width = inset.w, height = inset.h, vjust = 1)
# Fire
gg.america.burned.inset <- cowplot::ggdraw(gg.america.burned) +
  cowplot::draw_plot(gg.burned.fraction$`Pacific Northwest`, 
                     x = 0, y = 1, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.burned.fraction$California, 
                     x = 0, y = .6, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.burned.fraction$`Intermountain West`, 
                     x = .275, y = .4, 
                     width = inset.w, height = inset.h, vjust = 1) +
  cowplot::draw_plot(gg.burned.fraction$`Central America`, 
                     x = .775, y = .8, 
                     width = inset.w, height = inset.h, vjust = 1)

gg.asia.burned.inset <- cowplot::ggdraw(gg.asia.burned) +
  cowplot::draw_plot(gg.burned.fraction$Himalayas,
                     x = .15, y = 0.4,
                     width = inset.w, height= inset.h, vjust = 1) +
  cowplot::draw_plot(gg.burned.fraction$`Southeast Asia`,
                     x = .775, y = .935, 
                     width = inset.w, height = inset.h, vjust = 1)

# Precipitation
inset.w.p = .23
inset.h.p = .41
gg.america.precip.inset <- cowplot::ggdraw(gg.america.precip) +
  cowplot::draw_plot(gg.climatology$`Pacific Northwest`, 
                     x = 0, y = 1.01, 
                     width = inset.w.p, height = inset.h.p, vjust = 1) +
  cowplot::draw_plot(gg.climatology$California, 
                     x = .085, y = .64, 
                     width = inset.w.p, height = inset.h.p, vjust = 1) +
  cowplot::draw_plot(gg.climatology$`Intermountain West`, 
                     x = .245, y = .395, 
                     width = inset.w.p, height = inset.h.p, vjust = 1) +
  cowplot::draw_plot(gg.climatology$`Central America`, 
                     x = .71, y = .8, 
                     width = inset.w.p, height = inset.h.p, vjust = 1)

gg.asia.precip.inset <- cowplot::ggdraw(gg.asia.precip) +
  cowplot::draw_plot(gg.climatology$Himalayas,
                     x = .15, y = 0.4,  
                     width = inset.w.p, height= inset.h.p, vjust = 1) +
  cowplot::draw_plot(gg.climatology$`Southeast Asia`,
                      x = .7, y = .935, 
                     width = inset.w.p, height = inset.h.p, vjust = 1)
panel.labels <- tibble(
  label = as_factor(c('Location accuracy distribution\nfor burned and unburned sites by region',
                      'Locations of burned sites and\nregional distribution of burned fraction',
                      str_c('Landslide-triggering precipitation percentile\nand',
                            'regional seasonal precipitation (mm)')))
)
gg.panel.labels <- ggplot(panel.labels) +
  facet_wrap(~label, ncol = 1, switch = T) +
  theme(panel.background = element_blank(),
        strip.text = element_text(face="bold", size=7))

gg.panel.maps <- cowplot::plot_grid(
  cowplot::plot_grid(gg.america.region.inset, gg.asia.region.inset,
                   labels = c('a', 'b'), label_x = 1, hjust = 1, vjust = 1),
  legend.regions,
  cowplot::plot_grid(gg.america.burned.inset, gg.asia.burned.inset,
                     labels = c('c', 'd'), label_x = 1, hjust = 1, vjust = 1),
  legend.burned,
  cowplot::plot_grid(gg.america.precip.inset, gg.asia.precip.inset,
                     labels = c('e', 'f'), label_x = 1, hjust = 1, vjust = 1),
  legend.precip,
  ncol = 1, rel_heights = c(1, .3, 1, .2, 1, .3))
gg.panel <- cowplot::plot_grid(gg.panel.labels, gg.panel.maps, 
                               nrow = 1, rel_widths = c(.5, 6.5))
gg.panel

ggsave(filename = "../figures/panel_map.png", width=7, height=8)

Fire timing

This figure shows the timing of fires relative to landslides for all events that took place at a burned site.

calc.fire.timing <- function (modis.date, modis, slide) {
  modis.date %>%
    dplyr::filter(burn_days_before >= 0) %>%
    dplyr::group_by(OBJECTID) %>%
    dplyr::filter(burn_days_before == min(burn_days_before))  %>%
    right_join(modis %>% dplyr::filter(burned=='burned'), by='OBJECTID') %>%
    inner_join(slide.regions, by='OBJECTID') %>%
    left_join(slide %>% dplyr::select(OBJECTID, event_yday, season), by='OBJECTID') %>% 
    dplyr::mutate(burn_yday = yday(burn_date),
                  burn_season = factor(burn_season,
                                       levels = c('Spring', 'Summer', 'Fall', 'Winter'),
                                       ordered = T))
}

get.data('fire.timing', overwrite = F,
         calc.fire.timing, 
         modis.date, modis, slide)

Most landslides take place within one year of the fire. California stands out as having over a quarter of the landslides taking place more than a year after the fire.

fire.timing %>%
  dplyr::mutate(gt.1year = ifelse(burn_days_before > 2*365, 
                                  'gt.1year', 'lt.1year')) %>%
  dplyr::group_by(region, gt.1year) %>%
  dplyr::tally() %>%
  tidyr::spread(gt.1year, n, drop=F) %>%
  dplyr::mutate(prop.gt.1year = gt.1year / (gt.1year + lt.1year))
fire.to.landslide.x.breaks <- seq(season.breaks[4] - 4 * 365, season.breaks[3], 
                                  365/4)
fire.to.landslide.x.labels <- map(c(str_c('-', c(3, 2, 1), 'y '), ''), 
                                  ~c(str_c(., 'Winter'), 
                                     'Spring', 'Summer', 'Fall')) %>% as_vector

fire.facet.labels <- fire.timing %>%
  dplyr::group_by(region) %>%
  dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'], 
                                         event_yday - 365, event_yday)) %>%
  dplyr::ungroup() %>%
  dplyr::mutate(ymin = min(wrap_event_yday - burn_days_before),
                ymax = max(wrap_event_yday)) %>%
  dplyr::group_by(region) %>%
  dplyr::summarize(y = n() * 0.8, x = min(ymin) + (max(ymax) - min(ymin)) * 0.1) %>%
  dplyr::mutate(label = c('a', 'b', 'c', 'd', 'e', 'f', 'g'))

fire.timing %>%
  dplyr::group_by(region) %>%
  dplyr::mutate(wrap_event_yday = ifelse(event_yday > season.breaks['winter'], 
                                         event_yday - 365, event_yday)) %>%
  dplyr::arrange(wrap_event_yday - burn_days_before) %>%
  dplyr::mutate(i = row_number()) %>%
  ggplot() +
  facet_wrap(vars(region), scales = 'free_y', ncol = 2, shrink=F) +
  geom_segment(aes(x = wrap_event_yday - burn_days_before, xend = wrap_event_yday,
                   y = i, yend = i, 
                   color = burn_season)) +
  geom_point(aes(x = wrap_event_yday, y = i, shape='Landslide'), size=0.5) +
  geom_vline(data = tibble(x=seq(334-365*3, 365, 365)), 
             aes(xintercept = x, size = 'start_of_winter')) +
  geom_vline(aes(xintercept = 334 - 365, size = 'year_of')) +
  geom_rug(aes(y = i, color = burn_season)) +
  geom_rug(aes(x = event_yday - burn_days_before), sides = 't', length = unit(0.05, "npc")) +
  geom_label(data = fire.facet.labels,
             aes(x = x, y = y, label = label)) +
  theme_bw() +
  scale_x_continuous('Landslides and preceding fires relative to the year of the landslide',
                     minor_breaks = NULL,
                     breaks = as_vector(fire.to.landslide.x.breaks),
                     labels = fire.to.landslide.x.labels) +
  scale_color_manual('Fire season:',
                     values = c('Spring'='#19E667', 'Summer'='#E61998', 
                                'Fall'='#E6CD19', 'Winter'='#1932E6'),
                     guide = guide_legend(nrow=1, order = 2)) + 
  scale_y_continuous('Post-wildfire landslides ordered by fire timing relative to landslide', 
                     expand = c(0.1, 0.1)) +
  scale_size_manual('Relative timing:',
                    values = c('year_of' = 1, 'start_of_winter' = 0.5),
                    labels = c('year_of' = 'Start of year of landslide',
                               'start_of_winter' = 'Start of winter (December 1)'),
                    guide = guide_legend(nrow=2, order = 3)) +
  guides(shape = guide_legend(title = NULL, override.aes = c(size = 1.5), order = 1)) +
  theme(panel.grid.major.x = element_line(color = 'grey', size = 0.5),
        axis.text.x = element_text(hjust = 1, vjust = 1, angle = 90),
        axis.ticks = element_blank(),
        axis.text.y = element_blank(),
        legend.position = c(.75, .04),
        legend.direction = 'vertical',
        legend.box = 'vertical',
        legend.margin = unit(5, 'pt'),
        legend.box.margin = margin(5, 5, 5, 5),
        legend.box.background = element_rect(color = 'black', size = 2),
        legend.key.height = unit(10, 'pt') 
  )

ggsave('../figures/fire_to_landslide.png', width = 6.5, height = 8)

Analysis

CHIRPS vs Daymet

To validate the exclusion of zero-precipitation landslides, a comparison of CHIRPS and Daymet in California and Nevada only. Many of the zero-precipitation landslides measured by CHIRPS are not the same as the ones measured by Daymet, suggesting that these are measurement errors.

Approximately 14% of the landslides in California and Nevada were excluded on this basis.

join.chirps.daymet <- function (precip.chirps, precip.daymet) {
  precip.chirps %>%
    dplyr::inner_join(precip.daymet %>% dplyr::select(-days_before, -yday), 
                      by=c('OBJECTID', 'date', 'window'), 
                      suffix=c('.chirps', '.daymet'))
}

get.data('precip.chirps.daymet.ca.event', join.chirps.daymet, overwrite = F,
         precip.chirps.ca.event, 
         precip.daymet.ca.all)
get.data('precip.chirps.daymet.ca.all', join.chirps.daymet, overwrite = F,
         precip.chirps.ca.all, precip.daymet.ca.all)
precip.chirps.ca.event %>%
  filter(days_before==-1, window==7, !is.na(pctl)) %>%
  group_by(precip = ifelse(pctl > 0, 'non.zero.precip', 'zero.precip')) %>%
  tally() %>%
  spread(key = precip, value = n) %>%
  mutate(total = non.zero.precip + zero.precip,
         percent.zero.precip = zero.precip / total)
precip.chirps.daymet.ca.all %>%
  tidyr::drop_na(pctl.chirps, pctl.daymet) %>%
  ggplot() +
  geom_point(aes(x=pctl.chirps, y=pctl.daymet, color = 'all'), 
             alpha=0.2,shape = 16) +
  geom_point(data = precip.chirps.daymet.ca.event %>% 
               dplyr::filter(window == 7, days_before == -1 ) %>%
               dplyr::mutate(to.remove = ifelse(pctl.chirps==0, 'exclude', 'keep')),
             aes(x=pctl.chirps, y=pctl.daymet, color=to.remove)) +
  coord_equal() +
  theme_bw() +
  scale_color_manual('', 
                     values = c('keep'='black', 'exclude'='blue', 'all'='gray'),
                     labels = c('keep'='Included\nday-of-landslide\npreciptiation',
                                'exclude'='Excluded\nday-of-landlside\nprecipitation',
                                'all'='Precipitation\nfrom all dates')) +
  guides(color = guide_legend(keyheight = grid::unit(4, 'char'))) +
  labs(x='CHIRPS 7-day Precipitation Percentile',
       y='Daymet 7-day Precipitation Percentile')

ggsave('../figures/chirps_vs_daymet.png', width=6, height=4, dpi=600)

Precipitation percentile wilcox test

Comparison of antecedent precipitation time series among burned and unburned sites by region. Mann-Whitney tests are used to establish whether or not differences are statistically significant.

calc.burn.wilcox <- function (precip, pre=30, post=10, alternative='two.sided', 
                              source=NA, groupby=c('burned', 'days_before')) {
  precip %>%
    dplyr::filter(days_before <= pre, days_before >= -post) %>%
    dplyr::group_by_at(groupby) %>%
    calc.wilcox(alternative, source)
}
  
calc.wilcox <- function (precip, alternative='two.sided', source=NA) {
  wilcox <- precip %>%
    tidyr::nest() %>%
    spread(key = 'burned', value = 'data') %>%
    dplyr::mutate(wilcox = map2(burned, unburned, 
                         ~wilcox.test(.x$pctl, .y$pctl, 
                                      alternative=alternative))) %>%
    dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
    tidyr::unnest(wilcox.glance) %>%
    dplyr::mutate(signif.95 = p.value <= 0.05) %>%
    gather(key = 'burned', value = 'data', burned, unburned) %>%
    tidyr::unnest(data)
  if (!is.na(source)) wilcox <- wilcox %>% dplyr::mutate(label=source)
  wilcox
}

magnitude.regional.data <- function (slide.regions, precip, 
                                     pre = 30, post = 10, w = 7) {
  precip %>%
    dplyr::filter(days_before < pre, days_before > -post, window==w) %>%
    right_join(slide.regions, by = 'OBJECTID') %>%
    dplyr::left_join(modis, by = 'OBJECTID') %>%
    dplyr::filter(pctl > 0) %>%
    droplevels() %>%
    dplyr::select(region, burned, days_before, pctl) %>%
    calc.burn.wilcox(pre = pre, post = post, alternative = 'less',
                     groupby = c('region', 'burned', 'days_before'))
}

get.data('burn.wilcox.regions.event', overwrite = F,
         magnitude.regional.data,
         slide.regions,
         precip.chirps.ww.event)
facet.levels <- c('Rolling cumulative\nprecipitation percentile',
                  'Mann-Whitney p-value')

plot.burn.wilcox <- function (burn.wilcox.regions.event, slide.regions) {
  region.tally <- slide.regions %>% 
    left_join(modis) %>% 
    dplyr::group_by(burned=burn.cumulative > 0, region) %>% 
    tally
  
  burn.wilcox.regions.event %>%
    dplyr::mutate(signif.95 = ifelse(signif.95, 'significant', 'not.significant')) %>%
    ggplot() +
    facet_wrap(region~., ncol=1, strip.position = 'right') +
    geom_boxplot(aes(x = days_before, 
                     group=interaction(days_before, burned), 
                     y = pctl, 
                     fill = burned,
                     alpha = interaction(signif.95, burned)), 
                 outlier.shape=21)  +
    annotate("rect", xmin = -0.5, xmax = 0.5, ymin = 0, ymax = 1, alpha = .4) +
    geom_label(data=region.tally, 
               aes(x=-15, y=0.1+0.25*burned, 
                   label=str_c(ifelse(burned, 'Burned: ', 'Unburned: '), n)),
               hjust='right', label.size = 0) +
    geom_label(data = tibble(region = region.plot.order,
                             label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
                 dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
               aes(x = -15, y = 0.85, label = label), hjust = 'right') +
    scale_x_reverse(limits = c(31, -16), expand = expansion()) +
    scale_fill_manual('', 
                      values = c('burned' = '#d95f02', unburned = '#7570b3'),
                      labels = c('burned'='Burned', 'unburned'='Unburned'),
                      guide = guide_legend(direction = 'horizontal')) +
    scale_alpha_manual('Statistical significance at\n95% confidence',
                       values=c('not.significant.burned' = 0.4, 
                                'not.significant.unburned' = 0.4,
                                'significant.burned' = 1, 
                                'significant.unburned' = 1),
                       labels=c('not.significant.burned' = '', 
                                'not.significant.unburned' = 'Not Significant',
                                'significant.burned' = '', 
                                'significant.unburned' = 'Significant'),
                       breaks=c('not.significant.burned', 
                                'not.significant.unburned',
                                'significant.burned', 
                                'significant.unburned'),
                       guide=guide_legend(override.aes=list(fill=c('not.significant.burned' = '#d95f02', 
                                                                   'not.significant.unburned' = '#7570b3',
                                                                   'significant.burned' = '#d95f02', 
                                                                   'significant.unburned' = '#7570b3')),
                                          direction = 'horizontal')) +
    theme_bw() +
    labs(x='Days before the event', y='Rolling cumulative precipitation percentile') +
    theme(legend.position = 'bottom', legend.direction = 'vertical')
}

burn.wilcox.regions.event %>%
  plot.burn.wilcox(slide.regions)

ggsave('../figures/magnitude_regional.png', width = 11, height = 11)

Precipitation (mm) bootstrap sampling

Bootstrap samples are used here to develop baseline precipitation values against which to compare landslide-triggering precipitation. This method is an improvement over the original Wilcox tests because sample-size disparities between regions are reduced.

Bootstrap sampling

sample.bootstrap <- function(precip.subset, precip.all, n, samples, burned, region) {
  precip.all %>%
    dplyr::filter(burned==burned, region==region) %>%
    dplyr::mutate(precip.yday = yday(date)) %>%
    # dplyr::select precipitation from same location and date of year as subset
    dplyr::inner_join(precip.subset %>% 
                        dplyr::select(OBJECTID, event_yday), 
                      by=c('OBJECTID'='OBJECTID', 
                           'precip.yday'='event_yday')) %>%
    # Sample precipitation matching the site/day in bulk
    dplyr::group_by(OBJECTID) %>%
    tidyr::drop_na(pctl) %>%
    dplyr::sample_n(n, replace = T) %>%
    # Randomly split into individual samples
    dplyr::mutate(i = sample(1:samples, n, replace = T )) 
}

calc.burn.bootstrap <- function (precip.event, precip.all, 
                                 samples=100, per.subset=10000,
                                 pre=30, post=10, w=7) {
  groupby = c('region', 'burned', 'days_before')
  precip.all <- precip.all %>% 
    dplyr::filter(window==w) %>% 
    dplyr::select(-window)
  precip.event %>%
    rename(precip.date = date, precip.yday = yday) %>%
    dplyr::filter(days_before <= pre, days_before >= -post, window==w) %>%
    dplyr::select_at(c(groupby, 'OBJECTID', 'event_yday')) %>%
    # Tally each subset/site/day combination
    dplyr::group_by_at(groupby) %>%
    add_tally(name = 'sites.in.subset') %>%
    # Determine number to draws from each site that will 
    # roughly match sample sizes among subsets
    # Round up - number of draws must be an integer
    dplyr::mutate(n.per.site = ceiling(per.subset / sites.in.subset)) %>%
    # Take bootstrap samples for each subset/site/day combination
    dplyr::group_by_at(c(groupby, 'n.per.site', 'sites.in.subset')) %>%
    tidyr::nest() %>%
    dplyr::mutate(bootstrap.samples = pmap(list(data, n.per.site, 
                                                burned, region), 
                                           ~sample.bootstrap(..1, precip.all, 
                                                             ..2, samples, 
                                                             ..3, ..4)))
}
join.details <- function (precip, slide, slide.regions, modis) {
  precip %>% 
    dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday), 
                      by='OBJECTID') %>%
    dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
    dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned), 
                     by='OBJECTID')
}

get.bootstrap <- function (reg='All landslides', 
                           name = 'precip.chirps.bootstrap.all', 
                           ovw=F) {
  get.data(name, overwrite = ovw,
           calc.burn.bootstrap,
           precip.chirps.ww.event %>% 
             join.details(slide, slide.regions %>% filter(region == reg), modis),
           precip.chirps.ww.all %>% 
             join.details(slide, slide.regions %>% filter(region == reg), modis))
}

get.bootstrap('California', 'precip.chirps.bootstrap.cal', ovw = F)
get.bootstrap('Intermountain West', 'precip.chirps.bootstrap.imw', ovw = F)
get.bootstrap('Pacific Northwest', 'precip.chirps.bootstrap.pnw', ovw = F)
get.bootstrap('Central America', 'precip.chirps.bootstrap.cam', ovw = F)
get.bootstrap('Himalayas', 'precip.chirps.bootstrap.him', ovw = F)
get.bootstrap('Southeast Asia', 'precip.chirps.bootstrap.sea', ovw = F)
get.bootstrap(ovw = F)

precip.chirps.bootstrap <- precip.chirps.bootstrap.all %>%
  dplyr::bind_rows(precip.chirps.bootstrap.cal) %>%
  dplyr::bind_rows(precip.chirps.bootstrap.imw) %>%
  dplyr::bind_rows(precip.chirps.bootstrap.pnw) %>%
  dplyr::bind_rows(precip.chirps.bootstrap.cam) %>%
  dplyr::bind_rows(precip.chirps.bootstrap.him) %>%
  dplyr::bind_rows(precip.chirps.bootstrap.sea)

precip.chirps.bootstrap %>%
  dplyr::filter(days_before==0) %>%
  dplyr::select(-data) %>%
  dplyr::mutate(site.total = map_int(bootstrap.samples, ~nrow(ungroup(.))),
                n.per.site.real = site.total / sites.in.subset,
                avg.sample.size = map_dbl(bootstrap.samples, 
                                          ~group_by(., i) %>% 
                                            summarize(n = n()) %>% 
                                            pull(n) %>% mean)) %>%
  dplyr::select(-bootstrap.samples)

Summarize bootstrap samples

mann.whitney.p.value <- function (subset, bootstrap) {
  bootstrap %>%
    dplyr::group_by(i) %>%
    tidyr::nest() %>%
    # test if bootstrap samples are less than pre-landslide precipitation
    dplyr::mutate(wilcox = map(data,
                               ~wilcox.test(pull(., pctl), 
                                            pull(subset, pctl), 
                                            alternative='less'))) %>% 
    dplyr::mutate(wilcox.glance = map(wilcox, glance)) %>%
    dplyr::select(-data, -wilcox) %>%
    tidyr::unnest(wilcox.glance)
}
subset.precip <- function (precip, regions, modis, w=7, pre=30, post=10) {
  precip %>%
    dplyr::filter(window==7, days_before <= pre, days_before >= -post) %>%
    dplyr::select(OBJECTID, pctl, days_before) %>%
    dplyr::right_join(regions, by='OBJECTID') %>%
    dplyr::inner_join(modis, by='OBJECTID') %>%
    dplyr::group_by(region, burned, days_before) %>%
    tidyr::nest()
}

get_subset <- function (reg, bur, day, subsets) {
  row <- subsets %>%
    dplyr::filter(region==as.character(reg),
                  burned==bur,
                  days_before==day)
  row$data[[1]]
}

get.data('precip.chirps.ww.event.subsets', 
         subset.precip,
         precip.chirps.ww.event,
         slide.regions,
         modis)
summarize.bootstrap <- function (bootstrap.df, subsets, objective.function) {
  bootstrap.df %>%
    dplyr::mutate(objective = purrr::pmap(list(region, 
                                               burned, 
                                               days_before,
                                               bootstrap.samples), 
                                          ~objective.function(get_subset(..1, 
                                                                         ..2, 
                                                                         ..3, 
                                                                         subsets),
                                                              ..4))) %>%
    dplyr::select(-bootstrap.samples) %>%
    tidyr::unnest(objective) %>%
    dplyr::ungroup()
}
get.data('precip.chirps.bootstrap.wilcox', overwrite = F,
         summarize.bootstrap,
         precip.chirps.bootstrap %>% dplyr::select(-data),
         precip.chirps.ww.event.subsets,
         mann.whitney.p.value)
precip.chirps.bootstrap.wilcox %>% 
  dplyr::select(region, burned, days_before, p.value)

The bootstrap plot also shows the day-of-landslide kernel density estimates for all bootstrap samples against the landslide-triggering precipitation to illustrate the method. ### Plot bootstrap sampling

plot.burn.bootstrap <- function (bootstrap, var, name) {
  bootstrap %>%
    ggplot()  +
    facet_wrap(.~region, ncol=1, strip.position='right') +
    geom_vline(aes(xintercept=0)) +
    geom_boxplot(aes(x = days_before, 
                     group = interaction(days_before, burned), 
                     y = {{ var }}, 
                     fill = burned,
                     color = burned),
                 outlier.shape=21,
                 alpha = 0.6,
                 position = position_dodge(width=0.4, preserve = 'total'))  +
    annotate("rect", xmin = -11, xmax = -0.1, ymin = 0.05, ymax = .9999, 
             fill = 'white', alpha = 0.4) +
    annotate("rect", xmin = 0.1, xmax = 31, ymin = 0.05, ymax = .9999, 
             fill = 'white', alpha = 0.4)  +
    geom_hline(aes(yintercept=0.05), linetype='dashed') +
    geom_label(data = tibble(region = region.plot.order,
                             label = c('a', 'b', 'c', 'd', 'e', 'f', 'g')) %>%
                 dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
               aes(x = -11.5, y = 0.9999, label = label), 
               hjust = 'right', vjust = 'bottom')  +
    scale_x_reverse(expand = expansion(), limits = c(31, -12)) +
    scale_y_continuous(trans = probit_trans(),
                       breaks = bootstrap.y.breaks,
                       labels = bootstrap.y.labels) +
    scale_fill_burned +
    scale_color_burned +
    coord_trans(y = 'reverse', 
                ylim = c(10^-35, .9999)) +
    theme_bw() +
    labs(x='Days before the landslide', y=name) +
    theme(legend.position = 'bottom',
          panel.grid.minor.y = element_blank())
}

bootstrap.y.breaks <- c(10^-30,  10^-10,  0.05,   0.5,   0.9999)
bootstrap.y.labels <- c('1e-30', '1e-10', '0.05', '0.5', '1')
plot.bootstrap.wilcox <- precip.chirps.bootstrap.wilcox %>%
  plot.burn.bootstrap(p.value, 
                      'p-value of Wilcox tests between bootstrap samples and pre-landslide precipitation') 

plot.bootstrap.wilcox
plot.bootstrap.sample.kde <- precip.chirps.bootstrap %>%
  dplyr::filter(days_before==0) %>%
  dplyr::select(-data) %>%
  tidyr::unnest() %>%
  ggplot() +
  facet_grid(cols = vars(burned), rows = vars(region), scales = 'free_y') +
  geom_density(aes(x = pctl, group = i, color = burned, fill = burned), 
               alpha = 0.01) +
  geom_density(data = precip.chirps.ww.event %>% 
                 filter(days_before==0) %>% 
                 right_join(slide.regions, by = 'OBJECTID')%>%
                 inner_join(modis, by='OBJECTID'),
               aes(x = pctl, color = 'day_of', fill = 'day_of'), size = 1) +
  geom_label(data = cross_df(list(region = region.plot.order, 
                                  burned = c('burned', 'unburned'))) %>%
               dplyr::mutate(label = c('h', 'i', 'j', 'k', 'l', 'm', 'n', 
                                       'o', 'p', 'q', 'r', 's', 't', 'u')) %>%
               dplyr::mutate(region = factor(region, region.plot.order, ordered = T)),
             aes(x = .95, y = .05, label = label), 
             hjust = 'right', vjust = 'bottom') +
  scale_fill_manual('',
                    breaks = c('day_of', 'burned', 'unburned'),
                    values = c('day_of' = NA, burned.colors),
                    labels = c('day_of' = 'Day of landslide precipitation',
                               'burned' = 'Bootstrap samples, burned sites',
                               'unburned' = 'Bootstrap samples, unburned sites'),
                    guide = guide_legend(override.aes = list(alpha = 0.5))) +
  scale_color_manual('', 
                     breaks = c('day_of', 'burned', 'unburned'),
                     values = c('day_of' = 'black', burned.colors),
                     labels = c('day_of' = 'Day of landslide precipitation',
                                'burned' = 'Bootstrap samples, burned sites',
                                'unburned' = 'Bootstrap samples, unburned sites'),
                    guide = guide_legend(override.aes = list(alpha = 0.5))) +
  scale_x_continuous(breaks = seq(0, 1, 0.25),
                     labels = c('0', '0.25', '0.5', '0.75', '1')) +
  labs(y = 'Kernel density estimate', x = 'Precipitation percentile') +
  theme_bw() +
  theme(legend.position = 'bottom',
        legend.spacing.x = unit(15, "pt"),
        axis.text.y = element_blank(),
        axis.text.x = element_text(),
        axis.ticks.y = element_blank(),
        strip.text.x = element_blank(),
        panel.spacing.x = unit(15, 'pt')) +
  coord_cartesian(expand = F)
plot.bootstrap.sample.kde
plot.bootstrap.sample.kde.legend <- get_legend(plot.bootstrap.sample.kde)
plot.bootstrap.wilcox.legend <- get_legend(plot.bootstrap.wilcox)
cowplot::plot_grid(cowplot::plot_grid(plot.bootstrap.wilcox + 
                      theme(strip.text = element_blank(),
                            plot.margin = unit(c(5,5,5,5), 'pt'),
                            legend.position = 'none',
                            legend.direction = 'horizontal'), 
                    plot.bootstrap.sample.kde + 
                      theme(plot.margin = unit(c(5,5,5,2), 'pt'),
                            legend.position = 'none'),
                    rel_widths = c(11, 4)),
          cowplot::plot_grid(plot.bootstrap.wilcox.legend, 
                    plot.bootstrap.sample.kde.legend,
                    rel_widths = c(3, 7)),
          nrow = 2, rel_heights = c(10, 0.5))

ggsave('../figures/precip_bootstrap_wilcox+kde.png', width=13, height=11)

Seasonal density and precipitation frequency

Seasonal density

The kernel density estimate of the annual distribution of landslides is used to show shifts in seasonality relative to the day of the year.

slide.season <- slide %>%
  dplyr::select(OBJECTID, location_accuracy, season, event_yday) %>%
  dplyr::right_join(slide.regions, by='OBJECTID') %>%
  dplyr::left_join(modis, by='OBJECTID')

plot.seasonality <- bind_rows(slide.season %>% 
                                dplyr::mutate(event_yday = event_yday - 365),
                              slide.season,
                              slide.season %>% 
                                dplyr::mutate(event_yday = event_yday + 365)) %>%
  dplyr::group_by(region, burned) %>%
  tidyr::drop_na(burned) %>%
  tidyr::nest() %>%
  dplyr::mutate(kde = map(data, 
                          function (df) density(df$event_yday, 
                                                from=-100, to=450, 
                                                n = 512, bw = 21)),
                kde = map(kde, function (dens) tibble(event_yday = dens$x, 
                                                      density = dens$y))) %>%
  tidyr::unnest(kde) %>%
  dplyr::filter(event_yday >= 0, event_yday < 365) %>%
  ggplot() +
  facet_wrap(region~., ncol = 1, strip.position = 'right') +
  geom_area(aes(x = event_yday, y = density, group = burned, fill=burned),
            trim = T, alpha = 0.6, color = 'black', 
            position = position_identity()) +
  geom_text(data = tibble(text = c('h', 'i', 'j', 'k', 'l', 'm', 'n'),
                          hjust = c(-3.5, -10, -10, -4, -10, -2, -3.5),
                          region = region.plot.order,
                          x = -2.1*365, y = 0.50) %>%
               mutate(region = factor(region, region.plot.order, ordered = T)),
            aes(x = 0, y = .0035, label = text, hjust = hjust), size = 7) +
  coord_polar() +
  scale_fill_burned +
  scale_x_continuous(breaks = (season.breaks + 365/8) %% 365,
                     labels = c('W', 'Sp', 'Su', 'F'),
                     minor_breaks = season.breaks) +
  labs(y = '', x = ' ', fill = 'Season') +
  theme_bw() +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 0, vjust = 1),
        legend.position = 'bottom',
        panel.grid.minor.x = element_line(color = 'black'))

plot.seasonality

Precipitation frequency

Smoothed precipitation frequency is used to show shifts in the annual distribution of landslides relative to the annual precipitation seasonality.

n.years <- 2
calc.prop <- function (col) sum(!(col %in% 0)) / n()
calc.yday.frequency <- function (precip, slide, slide.regions, modis) {
  precip %>%
  dplyr::mutate(precip_yday = yday(date)) %>% 
  dplyr::select(OBJECTID, precip_yday, pctl) %>%
  dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_yday), 
                    by='OBJECTID') %>%
  dplyr::mutate(ydays_before = event_yday - precip_yday) %>%
  dplyr::select(OBJECTID, ydays_before, pctl) %>%
  dplyr::mutate(ydays_before = ydays_before %% 365) %>%
  dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
  dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned), 
                   by='OBJECTID') %>%
  dplyr::group_by(region, burned, ydays_before) %>%
  dplyr::summarize_at(vars(pctl), list(precip.yday.freq = calc.prop))
}
 
get.data('precip.chirps.yday.frequency',
         calc.yday.frequency,
         precip.chirps.ww.all,
         slide,
         slide.regions,
         modis)

calc.frequency <- function (precip, precip.yday.freq, 
                            slide, slide.regions, modis) {
  precip.chirps.ww.all %>%
    dplyr::transmute(OBJECTID, date = lubridate::as_date(date), pctl) %>%
    dplyr::inner_join(slide %>% dplyr::select(OBJECTID, event_date), 
                      by='OBJECTID') %>%
    dplyr::mutate(days_before = as.integer(event_date - date)) %>%
  dplyr::select(OBJECTID, days_before, pctl) %>%
    dplyr::filter(days_before < n.years * 365, 
                  days_before > -n.years * 365) %>%x
    dplyr::right_join(slide.regions, by = 'OBJECTID') %>%
    dplyr::left_join(modis %>% dplyr::select(OBJECTID, burned), 
                     by='OBJECTID') %>%
    dplyr::group_by(region, burned, days_before) %>%
    dplyr::summarize_at(vars(pctl), list(precip.day.freq = calc.prop)) %>%
    dplyr::mutate(ydays_before = days_before %% 365) %>%
    dplyr::left_join(precip.yday.freq, by=c('region', 'ydays_before', 'burned')) %>%
    dplyr::mutate(precip.freq.norm  = precip.day.freq / precip.yday.freq,
                  precip.freq.anomaly = precip.day.freq - precip.yday.freq,
                  precip.freq.lt.anomaly = precip.day.freq - mean(precip.day.freq)) %>%
    dplyr::ungroup()
}
get.data('precip.chirps.frequency', overwrite = F,
         calc.frequency,
         precip.chirps.ww.all,
         precip.chirps.yday.frequency,
         slide,
         slide.regions,
         modis)
plot.frequency <- function (data, column, y.lable, n.years, groupby=c('burned')) {
  data %>%
    dplyr::group_by_at(groupby) %>%
    dplyr::mutate(smooth = rollmean(!!as.name(column), 90, na.pad=TRUE)) %>%
    ggplot() +
    geom_line(aes(x = days_before, y = !!as.name(column), 
                  group = burned, color = burned, size = 'daily')) +
    geom_line(aes(x = days_before, y = smooth, group = burned, color = burned, 
                  size = 'rolling')) +
    geom_hline(data = data %>% 
                 dplyr::group_by_at(groupby) %>% 
                 summarize(mean = mean(!!as.name(column))),
               aes(yintercept = mean, color=burned), linetype = 'dashed') +
    geom_text(data = tibble(text = c('a', 'b', 'c', 'd', 'e', 'f', 'g'),
                            region = region.plot.order,
                            x = -1.95*365, y = 0.50),
              aes(x = x, y = y, label = text), size = 7,
              hjust = 'right') +
    scale_size_manual('',
                      values = c('daily' = 0.1, 'rolling' = 0.75),
                      labels = c('daily' = 'Daily frequency', 
                                 'rolling' = '90-day rolling average')) +
    scale_color_manual('', 
                      values = c('burned' = '#d95f02', unburned = '#7570b3'),
                      labels = c('burned'='Burned', 'unburned'='Unburned'),
                      guide = guide_legend(override.aes = list(size = 2))) +
    scale_x_reverse('Years before the landslide', 
                    breaks = seq(-365 * n.years, 365 * n.years, 365),
                    labels = function (breaks) breaks / 365,
                    expand = expansion()) +
    labs(y = y.lable) +
    theme_bw() +
    theme(legend.position = 'bottom')
}

plot.frequency.region <- precip.chirps.frequency %>%
  plot.frequency('precip.freq.lt.anomaly', 
                 'Precipitation frequency anomaly - long-term mean',
                 n.years=n.years, 
                 groupby=c('burned', 'region')) +
  facet_wrap(region~., ncol=1, strip.position = 'right')
plot.frequency.region
cowplot::plot_grid(
  cowplot::plot_grid(
    plot.frequency.region + 
      theme(strip.text.y = element_blank(),
            legend.position = 'none',
            plot.margin = unit(c(5, 1, 5, 5), 'pt')),
    plot.seasonality +
      theme(legend.position = 'none',
            plot.margin = unit(c(5, 5, 12, 0), 'pt')),
    rel_widths = c(10, 2.5)
  ),
  get_legend(plot.frequency.region),
  nrow = 2, rel_heights = c(10, 0.5)
)

ggsave('../figures/precip_frequency_seasonality.png', width=11, height=12)